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ABSTRACT 


This  report  provides  a  comprehensive  implementation  of  an  unstructured  mesh 
generation  method  that  refines  a  triangulated  grid  by  locally  bisecting  triangles 
on  their  longest  edge,  until  they  satisfy  a  given  local  condition.  The  method  is 
relatively  simple  to  implement,  has  the  capacity  to  quickly  generate  a  refined 
mesh  with  triangles  that  rapidly  change  size  over  a  short  distance,  and  does 
not  create  triangles  with  small  or  large  angles. 
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Mesh  Generation  via  Local  Bisection  Refinement  of 

Triangulated  Grids 

Executive  Summary 

A  mesh  is  a  collection  of  polygonal  (or  polyhedral)  elements  that  are  designed  to  approx¬ 
imate  a  geometric  domain.  Meshes  have  numerous  and  varied  applications  in  engineering 
and  scientific  computing,  such  as  the  numerical  solution  of  partial  differential  equations 
using  the  finite  element  method,  computer  aided  design,  rendering  computer  graphics,  and 
for  representing  terrain  surfaces  in  geographical  information  systems. 

This  report  emerged  from  a  study  of  path  planning  for  military  aircraft  conducting 
missions  in  hostile  environments,  where  the  objective  was  to  control  the  aircraft  such  that 
the  risk  posed  by  the  environment  was  minimised.  In  that  study,  mesh  elements  were 
smallest  in  regions  where  the  risk  to  the  aircraft  was  highest.  This  enabled  minimum  risk 
paths  to  be  accurately  captured  while  restricting  the  number  of  elements  in  the  mesh,  thus 
improving  the  efficiency  of  the  numerical  method  used  to  calculate  optimal  paths. 

This  report  documents  a  comprehensive  implementation  of  an  unstructured  mesh  gen¬ 
eration  method  that  refines  a  triangulated  grid  by  locally  bisecting  triangles  on  their 
longest  edge  until  they  satisfy  a  given  local  condition.  The  refined  mesh  consists  of  reg¬ 
ular  right-angled  triangles,  which  have  45,  45  and  90  degree  angles.  The  vertices  of  the 
triangles  are  ordered  such  that  minimal  computations  are  required  to  identify  their  longest 
edge  and  to  ensure  the  bisected  triangles  are  compatible,  that  is,  all  edges  are  shared  with 
at  most  one  other  triangle.  The  method  is  relatively  simple  to  implement,  has  the  capac¬ 
ity  to  quickly  generate  a  refined  mesh  with  triangles  that  rapidly  change  size  over  a  short 
distance,  and  does  not  create  triangles  with  small  or  large  angles. 

The  implementation  provided  by  this  report  has  the  following  desirable  features.  Mesh 
data  structures  are  chosen  to  enable  local  bisection  refinement  to  occur  in  constant  run¬ 
ning  time  and  with  minimal  computations.  A  local  refinement  condition  is  derived  that 
guarantees  the  local  triangulation  diameter  of  the  refined  mesh  will  obey  a  specified  bound; 
hence  the  numerical  error  of  computations  on  the  mesh  can  be  controlled  while  restricting 
the  number  of  triangles.  The  mesh  refinement  algorithm  is  tested  and  shown  to  achieve 
the  anticipated  linear  running  time  with  respect  to  the  number  of  triangles  in  the  refined 
mesh.  Employing  a  uniform  triangulated  grid  to  initialise  the  mesh  creates  virtual  buck¬ 
ets.  It  follows  that  point  location  can  be  accomplished  in  constant  running  time  without 
requiring  additional  memory. 


iii 


UNCLASSIFIED 


DSTO  -TR-3095 


UNCLASSIFIED 


THIS  PAGE  IS  INTENTIONALLY  BLANK 


IV 


UNCLASSIFIED 


UNCLASSIFIED 


DSTO-TR-3095 


Author 


Jason  Looker 

Joint  and  Operations  Analysis  Division 

Dr  Jason  Looker  joined  DSTO  in  2006  as  an  Operations  Re¬ 
search  Scientist  after  completing  a  BSc  (Honours)  and  PhD  in 
Mathematics  at  The  University  of  Melbourne.  He  has  under¬ 
taken  operations  analysis  in  support  of  Air  Force  Headquarters 
and  Project  AIR  9000  Phase  8  (Naval  Combat  Helicopter),  and 
is  leading  a  research  project  on  the  path  planning  of  military 
aircraft  through  hostile  environments. 


UNCLASSIFIED 


v 


DSTO  -TR-3095 


UNCLASSIFIED 


THIS  PAGE  IS  INTENTIONALLY  BLANK 


vi 


UNCLASSIFIED 


UNCLASSIFIED 


DSTO  -TR-3095 


Contents 

Notation  xi 

1  Introduction  1 

1.1  Local  Bisection  Refinement .  2 

1.2  Triangulation  Diameter .  4 

2  Mesh  Data  Structures  5 

2.1  Attributes .  5 

2.2  Memory  Allocation .  7 

3  Local  Refinement  Condition  9 

4  Coarse  Triangulation  10 

4.1  Construction .  10 

4.2  Mesh  Initialisation .  14 

5  Mesh  Refinement  15 

5.1  Bisecting  Triangles .  15 

5.2  Updating  Triangle  Neighbours  .  18 

5.3  Refinement  Algorithm .  20 

5.4  Computational  Performance .  21 

6  Adjacency  List  22 

7  Point  Location  23 

7.1  Initial  Triangle  Selection .  24 

7.2  Triangle  Containing  the  Query  Point .  25 

8  Conclusion  28 

References  29 


vii 


UNCLASSIFIED 


DSTO-TR-3095  UNCLASSIFIED 

Appendices 

A  Domains  other  than  [0, 1] 2  30 

B  Other  Cases  for  Bisecting  Triangles  and  Updating  their  Neighbours  33 


UNCLASSIFIED 


UNCLASSIFIED 


DSTO  -TR-3095 


Figures 

1  A  computational  mesh  generated  using  Maubach’s  method .  2 

2  A  coarse  triangulation  showing  the  assignment  of  indices  to  nodes  and  trian¬ 
gles,  and  the  grid  spacing  <5 .  11 

3  The  configuration  of  the  descendants  of  ti  for  the  case  k(tt)  =  1  and  O(L)  =  0. 

The  references  to  the  descendants  are  i  and  i,  and  a,  b ,  c  and  d  denote 
references  to  vertices.  Note  that  Nx  and  Nt  represent  their  states  immediately 
prior  to  bisecting  ti .  16 

4  The  configuration  of  triangle  neighbours  before  and  after  the  bisection  of 

compatibly  divisible  triangles  ti  and  tib  for  the  case  k{ti)  =  1,  0(U)  =  0  and 
0(tib)  =  0.  Here  a,  b,  c  and  d  represent  triangle  neighbours .  18 

5  The  scaled  average  running  time  of  the  implementation  of  Maubach’s  method 
documented  in  this  report,  versus  the  scaled  number  of  triangles;  a  line  of  best 

fit  is  shown  with  the  data  points .  22 

6  The  orientation  test:  does  d  lie  on,  to  the  left  of,  or  to  the  right  of  the  line 

defined  by  the  triangle  edge  (a,  &)?  The  arrows  indicate  the  directions  of  the 
vectors  used  to  perform  the  test.  The  test  is  applied  to  the  other  edges  in  the 
same  manner .  26 


UNCLASSIFIED 


IX 


DSTO  -TR-3095 


UNCLASSIFIED 


THIS  PAGE  IS  INTENTIONALLY  BLANK 


x 


UNCLASSIFIED 


UNCLASSIFIED 


DSTO  -TR-3095 


Notation 
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Zn 

Mn 

X 

1 1 X 1 1 

llxlll 

(ei,e2,...) 

{ei,e2,...} 

a  =  b  mod  i 

r-i 

L-J 

O(-) 
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Y{T) 

Nx 

Nt 
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h{x) 

h 

Mx) 

h{t) 

h{i) 

To 

8 

tb 

ib 


natural  numbers 

n-dimensional  space  of  integers 

n-dimensional  space  of  real  numbers 

a  generic  n-dimensional  vector,  x  =  (xi,  x2, . . . ,  xn) 

magnitude  (two-norm)  of  x,  |jx||  =  \Jy^a=i  X1 

one-nornr  of  x,  ||x||i  =  Y^=\  \xi\ 

an  ordered  list  with  elements  e,;  for  example,  a  vector 

an  unordered  list  with  elements  e*;  for  example,  a  set 

number  of  entries  in  a  data  structure 

modular  arithmetic 

ceiling  of  a  real  number 

floor  of  a  real  number 

“big  O”  Landau  symbol 

the  given  domain,  Slci2 

triangle  in  D 

triangulation  of  f 1 

set  of  nodes  that  comprise  T 

number  of  nodes  in  T 

number  of  triangles  in  T 

the  specified  maximum  grid  spacing  for  T 

triangulation  diameter 

local  triangulation  diameter  at  x  £  Y (T) 

the  specified  target  triangulation  diameter 

the  specified  target  local  triangulation  diameter  at  x  e  Y (T) 

the  specified  target  local  triangulation  diameter  at  t  £  T 

the  specified  target  local  triangulation  diameter  at  ti  G  T 

coarse  triangulation  of  D 

grid  spacing  of  To 

bisection  neighbour  of  t 

index  of  the  bisection  neighbour  of  U 

bisection  edge  length 

inverse  bisection  edge  length 
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level  of  refinement  of  t 

orientation  of  the  arrangement  of  the  vertices  of  t 
references  to  the  vertices  of  t,  V(t)  =  (vi,V2,vs) 
references  to  the  neighbours  of  t,  N(t )  =  (rai,  712,713) 
node  data  structure  (array) 
triangle  data  structure  (array) 
generic  data  structure  (array) 
ith  entry  in  V 
j th  entry  in  V[i] 

ith  to  jth  consecutive  entries  in  V 
j  th  to  Lth  consecutive  entries  in  V\i] 

v(u) 

N(U) 

0(u ) 

L(u) 

null  data  array  for  padding  J\f 
null  data  array  for  padding  T 
“true”  Boolean  data  type 
“false”  Boolean  data  type 
neighbourhood  of  x  £  Y (T) 

number  of  triangles  in  Tq  with  an  edge  on  any  one  side  of  [0,  l]2 
data  structure  for  references  to  triangle  vertices  in  T 
data  structure  for  references  to  triangle  vertices  in  To 
data  structure  for  references  to  triangle  neighbours  in  To 
reference  to  a  null  boundary  triangle,  (5  <  0 
determine  which  vertices  of  f *  form  the  bisection  edge 
stack  (LIFO  queue)  of  references  to  triangles  to  bisect 
adjacency  list  data  structure 
edges  of  f?; 

maximum  level  of  refinement  in  T 
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1  Introduction 

A  mesh  is  a  collection  of  polygonal  (or  polyhedral)  elements  that  are  designed  to  approx¬ 
imate  a  geometric  domain.  Meshes  have  numerous  and  varied  applications  in  engineering 
and  scientific  computing,  such  as  the  numerical  solution  of  partial  differential  equations 
using  the  finite  element  method,  computer  aided  design,  rendering  computer  graphics,  and 
for  representing  terrain  surfaces  in  geographical  information  systems. 

This  report  emerged  from  a  study  of  path  planning  for  military  aircraft  conducting 
missions  in  hostile  environments,  where  the  objective  was  to  control  the  aircraft  such  that 
the  risk  posed  by  the  environment  was  minimised.  In  that  study,  mesh  elements  were 
smallest  in  regions  where  the  risk  to  the  aircraft  was  highest.  This  enabled  minimum  risk 
paths  to  be  accurately  captured  while  restricting  the  number  of  elements  in  the  mesh, 
thus  improving  the  efficiency  of  the  numerical  method  used  to  calculate  optimal  paths 
[Looker  2013]. 

Meshes  are  either  structured  or  unstructured.  A  node  (element  vertex)  of  a  structured 
mesh  is  referenced  by  a  tuple  of  indices  (one  for  each  spatial  dimension),  and  adjacent  nodes 
are  found  by  translating  each  index  of  the  tuple  by  unity.  A  node  of  an  unstructured  mesh 
is  referenced  by  a  single  index,  and  topological  attributes  of  the  mesh  are  stored  in  separate 
data  structures  to  enable  adjacent  nodes  and  elements  to  be  found  [Hjelle  &  Dashlen  2006]. 
Note  that  the  characterisation  of  a  mesh  as  either  structured  or  unstructured  relates  to 
how  the  mesh  is  stored  in  memory;  it  does  not  relate  to  the  geometric  structure  of  the 
mesh.  However,  structured  meshes  must  have  an  underlying  geometric  structure,  whereas 
this  may  or  may  not  be  the  case  for  unstructured  meshes  [Nelson  2014]. 

Structured  meshes  are  simpler  to  implement  and  faster  to  generate  compared  with 
unstructured  meshes.  However,  unstructured  meshes  allow  greater  flexibility  in  the  choice 
of  element,  and  can  attain  a  rapid  change  in  element  size  over  a  shorter  distance  using 
fewer  elements  than  a  structured  mesh,  resulting  in  more  efficient  computations  on  the 
mesh  [Shewchuk  19976,  Shewchuk  2012,  Nelson  2014], 

This  report  documents  a  comprehensive  implementation  of  the  unstructured  mesh 
generation  method  of  Maubach  [1995],  focussing  on  the  case  of  a  two-dimensional  mesh; 
although  the  method  can  be  applied  in  higher  dimensions  [Maubach  1995,  Arnold,  Mukher- 
jee  &  Pouly  2000].  Maubach’s  method  refines  a  triangulation  (a  mesh  with  triangular 
elements)  by  locally  bisecting  triangles  on  their  longest  edge  until  they  satisfy  a  given 
local  condition.  The  refined  mesh  consists  of  regular  right-angled  triangles,  which  have 
45,  45  and  90  degree  angles,  as  shown  in  Figure  1.  The  vertices  of  the  triangles  are  ordered 
to  allow  their  longest  edge  to  be  identified  without  any  computations1  and  to  ensure  the 
bisected  triangles  are  compatible,  that  is,  all  edges  are  shared  with  at  most  one  other 
triangle. 

Maubach’s  method  is  relatively  simple  to  implement,  has  the  capacity  to  quickly  gen¬ 
erate  a  refined  mesh  with  triangles  that  rapidly  change  size  over  a  short  distance,  and 
does  not  create  triangles  with  small  or  large  angles  (which  may  cause  numerical  difficul¬ 
ties).  However,  all  refined  elements  are  regular  right-angled  triangles  and  this  may  be  too 

lrTlie  phrase  “without  any  computations”  is  to  be  interpreted  as  “without  any  searching,  sorting  or 
numerical  computations”  throughout  this  report. 
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Figure  1:  A  computational  mesh  generated  using  Maubach’s  method. 


restrictive  for  some  applications  [Shewchuk  2012].  Furthermore,  other  techniques  such  as 
Delaunay  mesh  generation  [Cheng,  Dey  &  Shewchuk  2012]  produce  size-optimal  meshes 
that  have  fewer  triangles  without  losing  numerical  precision,  and  hence  more  efficient 
computations  can  be  performed  on  the  resulting  mesh. 

The  introduction  continues  with  a  review  of  local  bisection  refinement  where  the  key 
algorithms  of  Maubach’s  method  are  presented,  and  concludes  with  a  brief  discussion  of  the 
triangulation  diameter.  Section  2  examines  the  data  structures  chosen  to  store  the  mesh. 
A  condition  for  local  mesh  refinement  is  derived  in  Section  3  with  the  aim  of  controlling 
the  numerical  error  of  computations  performed  on  the  mesh.  Initial  mesh  construction 
for  the  case  of  a  uniform  triangulated  grid  is  presented  in  Section  4.  Section  5  contains 
the  implementation  of  Maubach’s  method  that  is  the  principal  subject  of  this  report.  An 
algorithm  for  constructing  an  adjacency  list  of  the  nodes  from  the  refined  mesh  is  presented 
in  Section  6.  Determining  a  triangle  that  contains  a  given  point  is  a  fundamental  operation 
on  a  mesh,  and  is  known  as  point  location;  this  is  discussed  in  Section  7. 


1.1  Local  Bisection  Refinement 

In  this  section,  some  elementary  notation  and  terminology  is  introduced,  followed  by  a 
review  of  local  bisection  refinement. 
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Let  H  C  K2  be  a  given  domain.  A  triangle  in  12  is  denoted  by  t  where  x,;  G  12  are  the 
vertices  of  t.  A  triangulation  T  of  12  is  defined  by2 

T  =  {t\  t  do  not  overlap,  are  compatible,  and  all  x,;  G  11} . 

The  nodes  of  T  are  defined  to  be  the  set  of  vertices  that  comprise  all  t  G  T.  A  more  precise 
definition  of  T  can  be  found  in  Hjelle  &  Dashlen  [2006]. 

Maubach’s  method  begins  with  a  coarse  triangulation  To  of  12,  where  all  to  G  To  are 
identical.  The  number  of  bisections  of  to  required  to  generate  a  descendant  t  G  T  is  the 
level  of  refinement  of  t,  given  by  T(t);  the  level  of  refinement  of  all  to  is  defined  to  be  zero. 
Maubach’s  method  generates  a  refined  T  by  locally  bisecting  triangles  on  their  longest 
edge,  which  shall  be  referred  to  as  the  bisection  edge. 

Local  bisection  refinement  is  accomplished  by  two  procedures:  the  first  bisects  a  trian¬ 
gle  and  maintains  an  ordering  of  the  descendant’s  vertices  to  enable  their  bisection  edges 
to  be  subsequently  identified;  and  the  second  procedure  guarantees  the  descendants  are 
compatible. 

Initially  all  t  G  To  have  their  vertices  arranged  anticlockwise  such  that  the  first  and 
third  vertices  form  the  bisection  edge.  BisectTriangle  (Algorithm  1)  bisects  a  given 
triangle  t  and  creates  descendants  t\  and  t2 ■  The  order  of  the  vertices  of  t\  and  f2  is  chosen 
to  ensure  their  bisection  edges  can  be  immediately  identified,  without  any  computations, 
if  they  are  subsequently  bisected.  However,  BisectTriangle  is  not  sufficient  to  generate 
a  refined  T  as  mesh  incompatibilities  will  be  introduced. 


Algorithm  1:  BisectTriangle 


Input 

:  t  G  T  where  t  has  vertices  (xi,X2,X3) 

Result:  descendants  t\  and  2 2  are  created 

1 

k 

«-  2 

—  (T(2)  mod  2) 

2 

z 

1 

2 

(Xl  +  Xfc+i) 

3 

if  k  = 

1  then 

4 

t\  - 

(xi,z,x3) 

5 

t2  - 

(x2,z,x3) 

6 

else 

7 

t\  - 

(xi,x2,z) 

8 

t2  - 

(x2,x3,z) 

9 

L(h)  < 

—  L(t)  +  1 

10 

—  L(t )  +  1 

A  record  of  triangle  neighbours  is  maintained  to  facilitate  the  avoidance  of  mesh  in¬ 
compatibilities.  Two  triangles  are  neighbours  if  they  share  an  edge.  Compatibly  divisible 
triangles  are  neighbours  that  share  their  bisection  edge.  The  neighbour  of  t  on  its  bisection 
edge  is  called  the  bisection  neighbour  and  denoted  by  Incompatibilities  are  avoided 
if  two  compatibly  divisible  triangles  are  bisected  simultaneously.  If  L(t)  =  L(tb)  then  2 

2  The  symbol  T  may  denote  a  generic  triangulation  or  the  state  of  a  triangulation  being  refined  using 
Maubach’s  method. 


UNCLASSIFIED 


3 


DSTO  -TR-3095 


UNCLASSIFIED 


and  4  are  compatibly  divisible,  otherwise  t  will  be  compatibly  divisible  with  one  of  the 
descendants  of  tb  [Maubach  1995,  Theorem  5.1]. 

This  is  exploited  in  the  recursive  algorithm  Ref  ineTriangle  (Algorithm  2)  to  com¬ 
patibly  refine  a  given  triangle;  the  recursion  depth  of  Ref  ineTriangle  is  bounded  by  the 
maximum  level  of  refinement  in  T  [Maubach  1995].  Ref  ineTriangle  calls  itself  repeatedly 
on  a  sequence  of  triangles  until  a  compatibly  divisible  triangle  is  found.  This  sequence  of 
triangles  is  then  bisected  in  reverse  order  to  preserve  compatibility. 


Algorithm  2:  Ref  ineTriangle 


Input:  i  £  T 

Result:  t  is  compatibly  refined 


1 

2 

3 

4 

5 

6 

7 

8 


tb  bisection  neighbour  of  t 
if  L(t )  =  L(tb)  then 
BisectTriangle(f) 

BisectTriangle(4) 

foreach  compatibly  divisible  triangle  t!  visited  by  Ref  ineTriangle  do 
t'b  <—  bisection  neighbour  of  t! 

BisectTri  angle  (0) 

BisectTri  angle  (t'b) 


9  else 

10  add  t  to  the  sequence  of  triangles  visited  by  Ref  ineTriangle 

11  Ref ineTriangle(tfe) 


Maubach’s  method  has  linear  running  time  with  respect  to  the  number  of  trian¬ 
gles  in  the  final  refined  mesh.  This  follows  from  the  bound  on  the  recursion  depth  of 
Ref  ineTriangle  that  only  depends  on  the  maximum  level  of  refinement  in  T,  and  since 
Maubach’s  method  is  based  on  bisection  refinement. 


1.2  Triangulation  Diameter 

A  key  attribute  of  triangulations  relevant  to  scientific  computing  applications  is  the  trian¬ 
gulation  diameter  h,  as  the  numerical  error  of  computations  performed  on  triangulations 
can  often  be  estimated  in  terms  of  h.  The  triangulation  diameter  of  T  is  defined  to  be  the 
maximum  distance  between  adjacent  nodes  in  the  mesh;  two  nodes  are  adjacent  if  they 
form  an  edge  of  a  t  €  T. 

In  this  report  the  coarse  triangulation  To  is  defined  to  be  a  uniform  triangulated  grid, 
where  all  t  £  To  are  identical  regular  right-angled  triangles  with  vertices  coinciding  with 
the  nodes  of  a  square  grid  with  spacing  5.  The  bisection  edge  length  lb  is  given  by 

4(5,0  =  V2{1~1)8,  (1) 

for  t  6  T  with  l  =  L(t).  Equation  (1)  follows  from  the  Pythagorean  theorem  and  the 
definition  of  bisection  refinement.  It  follows  that  the  triangulation  diameter  of  T  generated 
using  Maubach’s  method  is  given  by  h  =  s/2  5,  since  L(t)  =  0  for  all  t  €  T> 
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2  Mesh  Data  Structures 

Maubach’s  method  generates  an  unstructured  mesh  and  hence  data  structures  must  be 
defined  with  particular  fields  to  store  the  geometric  and  topological  attributes  of  the  mesh. 
Many  types  of  mesh  data  structure  exist:  some  minimise  the  memory  required  to  store  the 
mesh;  others  minimise  the  speed  of  performing  certain  operations  on  the  mesh.  Ultimately, 
the  choice  of  a  particular  mesh  data  structure  is  highly  context  dependent.  Refer  to  Hjelle 
&  Daehlen  [2006]  for  a  general  discussion  of  data  structures  for  triangulations. 

It  can  be  seen  from  Algorithms  1  and  2  that  the  following  mesh  data  must  be  stored: 

•  nodes  (which  form  the  triangle  vertices) 

•  references  to  triangle  vertices 

•  references  to  triangle  neighbours 

•  refinement  levels. 

The  orientation  of  the  arrangement  of  triangle  vertices  is  also  required  to  be  stored  for 
updating  references  to  triangle  neighbours,  the  discussion  of  which  is  postponed  until 
Section  5.  The  orientation  data  is  also  used  for  point  location,  which  is  discussed  in 
Section  7. 

For  each  call  to  Ref  ineTriangle,  node  and  triangle  data  must  be  found,  new  nodes 
and  triangles  created,  and  some  existing  triangles  must  have  their  data  updated.  Therefore 
three  operations  are  required  to  be  performed  on  the  mesh  data  structures:  find  entries; 
add  new  entries;  and  update  existing  entries. 

A  natural  choice  of  mesh  data  structure  is  an  array,  as  node  and  triangle  data  can  be 
referenced  using  their  unique  array  indices,  and  looking  up  a  particular  entry  in  an  array 
is  very  fast.  However,  searching  and  adding  entries  to  an  array  are  expensive  operations, 
especially  if  the  array  is  very  large.  Also  note  that  rectangular  arrays  with  entries  of  the 
same  type  (integers,  reals,  . . . )  are  more  computationally  efficient  to  store  and  access, 
compared  with  non- rectangular  (ragged)  arrays  with  entries  of  mixed  type. 

These  observations  motivate  the  choice  of  mesh  data  structures  employed  in  the  im¬ 
plementation  of  Maubach’s  method  discussed  in  this  report. 


2.1  Attributes 

Let  Af  be  the  data  structure  for  nodes,  where  Af  is  an  array  with  entries  A f[i\  =  x,;  £  D. 
The  vertices  of  t  6  T  are  denoted  by  V (t)  and  specified  by  references  to  entries  in  Af;  for 
example,  if  t  has  the  vertices  (X95,  X154,  X3)  then  V(t)  =  (95,154,3).  For  the  general  case 
let  V(t)  =  (vi,  V2 ,  U3 )  where  Vi  £  N  are  references  to  entries  in  A f.  The  Vi  are  ordered 
according  to  Maubach’s  method;  refer  to  Algorithm  1. 

Let  T  be  the  data  structure  for  triangles,  where  T  is  an  array  with  entries  T[i\  £  Z8 
that  contain  data  associated  with  the  ith  triangle  in  T.  The  triangle  neighbours  of  t  £  T 
are  denoted  by  N(t)  =  (ni,  ri2,  ny)  where  nt  £  Z.  If  m  >  0  then  n*  references  a  triangle  in 
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T.  otherwise  if  rii  <  0  then  n*  references  a  null  boundary  triangle.3  The  n*  are  ordered  to 
facilitate  the  identification  of  neighbours  without  any  computations: 


•  ni  references  the  bisection  neighbour 

•  neighbours  are  arranged  anticlockwise  about  t 

•  if  rii  <  0  then  t  has  an  edge  on  the  boundary  of  f I  that  would  be  shared  with  m,  if 
it  was  not  a  reference  to  a  null  boundary  triangle. 


Let  0(t )  E  {0, 1}  represent  the  orientation  of  the  arrangement  of  the  vertices  of  t  E  T, 
where  0(t )  =  0  indicates  the  vertices  are  arranged  clockwise  and  0(t)  =  1  indicates  the 
vertices  have  an  anticlockwise  arrangement.  Further  discussion  of  0{t )  is  postponed  until 
Section  5. 

The  zth  entry  of  T  is  given  by 

'I'i'i]  —  (vil )  ^ilj  Oi,  li) , 

where  o*  =  ()(tt)  and  li  =  L{ti).  Let  T[i,j]  refer  to  the  jth  entry  in  T[i\  for  1  ^  j  ^  8, 
and  T[i,  j  :  k]  refer  to  the  jth  to  kth  consecutive  entries  of  T[i\  for  1  ^  j  <  k  ^  8.  It 
follows  that 


T[i,l:3\  =  V(U), 
T[i,4  :  6]  =  N(ti), 
T[i,7]  =  0(U), 
T[i,8\  =  L(ti). 


The  attributes  of  M  and  T  facilitate  a  computationally  efficient  implementation  of 
Maubach’s  method: 


•  M  and  T  are  rectangular  arrays  each  with  entries  of  the  same  type 

•  the  order  of  the  entries  in  T[i]  enables  references  to  triangle  vertices  and  neighbours 
to  be  identified  without  any  computations 

•  the  compact  structure  of  T[i\  minimises  the  need  to  append  triangle  data  to  multiple 
data  structures. 

In  fact,  it  will  be  shown  in  Section  2.2  that  appending  data  to  M  and  T  can  almost  be 
eliminated. 

3The  case  rn  =  0  will  be  discussed  in  Section  2.2. 
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2.2  Memory  Allocation 

The  triangulation  T  is  initialised  to  To  and  stored  in  J\f  and  T,  then  calls  to  Ref  ineTriangle 
locally  refine  T  resulting  in  new  nodes  and  triangles  being  created.  One  option  is  to  ap¬ 
pend  new  node  and  triangle  data  to  N  and  T ;  however  this  is  an  expensive  operation.  A 
more  efficient  alternative  is  to  estimate  the  sizes  of  M  and  T  necessary  to  store  the  final 
refined  T,  and  pad  out  N  and  T  with  the  required  amount  of  null  data  after  initialisation. 
These  estimates  should  be  upper  bounds  as  this  will  eliminate  the  requirement  to  append 
data  to  M  and  T ;  the  excess  padding  can  be  removed  once  the  refinement  of  T  is  complete. 

Maubach’s  method  refines  T  according  to  a  given  local  condition,  which  in  this  report 
is  assumed  to  depend  on  a  function  that  specifies  the  target  bisection  edge  length  bt  for 
each  t  £  T.  To  estimate  the  number  of  triangles  in  the  final  refined  T,  L(t)  will  be 
estimated  for  each  t  £  To  using  bt-  Then  YlteTu  2 L ^  triangles  would  be  created  if  each 
t  £  To  and  all  their  descendants  were  bisected  L(t)  times.  The  inverse  bisection  edge 
length  £^(8,  bt)  is  defined  to  return  the  minimum  L(t )  such  that  £b(5,L(t))  ^  bt  and  is 
given  by5 

4_1(W  =  [21og2(72  6/bt)], 
which  follows  from  Equation  (1). 

The  number  of  triangles  in  the  final  refined  T  is  denoted  by  Nf  and  can  be  estimated 
in  terms  of  £^1: 

Nf.  «  ^  2^~1('5’6‘).  (2) 

teTo 

It  can  be  shown  that  the  number  of  nodes  and  triangles  in  a  triangulation,  and  hence  the 
number  of  entries  in  J\f  and  T,  are  related  via  [Hjelle  &  Daehlen  2006,  Lemma  1.1] 

|T|  ~  2  |AT|  as  |jV|  — ►  oo.  (3) 

Hence  the  number  of  nodes  in  the  final  refined  T  will  be  estimated  using 

\nI  ,  (4) 

with  N't  given  by  Equation  (2). 

Since  A f[i\  £  M2  and  T[i]  £  Z8,  the  arrays  of  null  data  used  to  pad  out  M  and  T  are 
defined  to  be 

VM  =  ((0.0,  0.0),...), 

Vr  =  ((0,0, 0,0, 0,0, 0,0),...). 

At  initialisation  of  Maubach’s  method,  T  <—  To,  the  mesh  data  is  stored  in  J\f  and  T,  and 
then 

A  <—  (AT[l] , AT[2] , . . .  ,AA[NX],TV[1],  •  •  •  ,TVAx  ~  ^x]),  (5) 

T  e-  (T[l] ,  T[2] , . . . ,  T[Nt } ,  Vr  [1],  •  • . ,  Vr  [Ntf  -  (6) 

4The  refinement  condition  is  discussed  in  Section  3. 

5 Note  that  £ is  not  a  true  mathematical  inverse. 
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where  Atf  and  A£  are  given  by  Equations  (2)  and  (4),  respectively,  and  Ax  is  the  number 
of  nodes  and  Nt  the  number  of  triangles  in  T. 

Equations  (2)  and  (4)  represent  informal  upper  bounds  on  A |  and  Ah,  and  therefore 
J\f  and  T  as  defined  by  Equations  (5)  and  (6)  may  have  insufficient  padding  to  store  the 
final  refined  T.  Consequently  the  residual  padding  must  be  checked  before  each  iteration 
of  Maubach’s  method.  Ref  ineTriangle  is  a  recursive  algorithm  and  hence  the  number 
of  nodes  and  triangles  created  per  call  cannot  be  known  in  advance,  and  for  this  reason 
the  residual  padding  must  be  checked  before  bisection  is  performed  by  BisectTriangle. 
Observe  from  Algorithm  2  that  BisectTriangle  only  bisects  compatibly  divisible  triangles 
and  their  bisection  neighbours,  creating  one  node  and  two  triangles  per  call.6  Therefore  if 

W\>NX  +  1,  (7) 

and 

|T|  >Nt  +  2,  (8) 

then  there  is  sufficient  padding  in  J\f  and  T  for  at  least  one  call  to  BisectTriangle.  In 
fact,  it  is  now  shown  that  Equation  (7)  implies  Equation  (8). 

Let  Mq-  ^  2  be  a  given  number  of  additional  null  triangles  that,  if  required,  can  be 
used  to  pad  out  T.  Motivated  by  Equation  (4),  \_Mq-/2\  null  nodes  will  be  used  to  pad 
out  M  as  required.  Accordingly,  \Pr\  ^  2  \Vj\f\  for  all  iterations  of  Maubach’s  method.  It 
follows  that  if  Equation  (7)  is  true,  then  for  all  iterations  of  Maubach’s  method, 

JVt  +  2<|T|-2|7V|  +  2, 

=  \T\  —  2  (|A7|  —  A^)  +  2, 

=  \T\  —  2  |A7|  +  2  (Ax  +  1), 


<i  n 

and  therefore  Equation  (8)  is  also  true. 

The  preceding  discussion  motivates  the  definition  of  CheckMemory  presented  in  Algo¬ 
rithm  3.  Calling  CheckMemory  before  bisection  is  performed  by  BisectTriangle  guaran¬ 
tees  there  will  be  sufficient  null  data  in  M  and  T  to  store  new  nodes  and  triangles. 


Algorithm  3:  CheckMemory 

Result:  A7  and  T  contain  enough  padding  to  call  BisectTriangle  on  a  triangle 
and  its  compatibly  divisible  neighbour 

1  if  |A7|  <  Ax  +  1  then 

2  M  <-  (A7[l], A7[2], . . .  ,A7[AX],7V[1],. . .  ,VM[  [MT/2\  ]) 

3  T  <-  (T[1],T[2],  ...,T[Nt],Vr[l],..  ■  ,VT[MT}) 


6  Bisecting  a  triangle  that  has  its  bisection  edge  on  the  boundary  of  creates  one  node  and  one  triangle 
per  call  to  BisectTriangle;  this  issue  will  be  considered  in  Section  5. 
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3  Local  Refinement  Condition 

The  aim  of  this  section  is  to  derive  a  local  refinement  condition  that  returns  True  if 
Ref  ineTriangle(f)  must  be  called,  or  False  otherwise;  Maubach’s  method  will  terminate 
when  this  condition  returns  False  for  all  t  £E  T.  This  local  refinement  condition  will 
be  derived  such  that  the  distance  between  neighbouring  nodes  in  the  final  refined  T  will 
be  bounded  by  a  specified  target,  with  the  goal  of  controlling  the  numerical  error  of 
computations  performed  on  T,  while  curtailing  Nf. 

First,  the  local  triangulation  diameter  is  defined.  The  neighbourhood  ./F(x)  of  a  node 
x  in  T  is  given  by 

o/K(x)  =  {y  6  O  |  there  exists  a  i  £  T  such  that  x£i  and  y  G  t},  (9) 

where  the  notation  z  £E  t  is  to  be  interpreted  as  “z  is  a  vertex  of  t” .  The  local  triangulation 
diameter  /t(x)  is  defined  in  terms  of  ,yU(x): 

/i(x)  =  max  1 1 x  —  y  1 1  •  (10) 

ye^F(x) 

Let  y(T)  represent  the  set  of  nodes  that  comprise  T,  that  is, 

y  (T)  =  {x  £  9  |  there  exists  a  t  G  T  such  that  x£f}. 

Equation  (10)  provides  a  more  general  definition  of  the  triangulation  diameter  h  than  that 
given  in  Section  1.2: 

h  =  max  /i(x).  (11) 

xer(T) 

Let  h(x)  be  the  specified  target  local  triangulation  diameter  for  all  x  G  y(T).  Since  the 
local  refinement  condition  will  be  defined  in  terms  of  triangles,  define  h(t)  via 

h(t)  =  minlt(x),  (12) 


for  all  i  £  T.  Note  that  h(t)  is  equivalent  to  bt  discussed  in  Section  2.2. 

The  local  refinement  condition  is  given  by 

4  (6,L(t))>h(t),  (13) 

where  5  is  the  grid  spacing  of  To,  and  recall  that  Maubach’s  method  will  terminate  when 
Equation  (13)  returns  False  for  all  t  £  T.  Let  T  denote  a  triangulation  generated  using 
Maubach’s  method.  Equation  (13)  ensures  T  will  obey 

h{x)  4  h(x)  for  all  x  €  T(T).  (14) 

To  verify  Equation  (14),  first  let 


^(x)  =  {t  e  T  |  x£f}, 
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and  choose  any  x  €  Y(T).  Then 

h(x)  =  max  1 1  x  —  yll. 

V  ye^(x) 


^  max  £b(5,  L(t))  (property  of  local  bisection  refinement), 
te^x) 


^  max  h(t)  (by  the  negation  of  Equation  (13)), 

te^(x) 


=  max  min/i(y)  (by  Equation  (12)), 

te^x)  ye* 


^  /i(x)  (as  xU  for  all  t  €  7^(x)). 

The  given  target  triangulation  diameter  h  is  analogously  defined  to  h  in  Equation  (11). 
Therefore  Equation  (13)  enables  h  of  the  final  refined  T  to  be  controlled  using  h,  since 
Equation  (14)  implies  h  ^  h,  thereby  controlling  the  numerical  error  of  computations 
performed  on  T. 


4  Coarse  Triangulation 


Maubach’s  method  begins  with  a  coarse  triangulation  To  of  Q.  To  simplify  the  discussion 
fl  is  defined  to  be  [0, 1] 2;  non-square  domains  and  domains  with  holes  are  considered  in 
Appendix  A.  The  aim  of  this  section  is  to  present  algorithms  for  constructing  Tq. 


Maubach’s  method  can  simply  be  initialised  using  a  To  containing  two  triangles  that 
are  compatibly  divisible.  However  in  this  report  To  is  chosen  to  be  a  uniform  triangulated 
grid,  where  all  t  £  To  are  identical  regular  right-angled  triangles  with  vertices  coinciding 
with  the  nodes  of  a  square  grid  with  spacing  5.  In  Section  7  this  choice  will  be  shown  to 
facilitate  efficient  point  location.  The  assignment  of  indices  to  nodes  and  triangles  in  To 
is  shown  in  Figure  2. 


The  specified  maximum  grid  spacing  5  for  T  is  related  to  5  via 

i-i 


S  = 


6  1 


(15) 


This  definition  of  5  implies  that  h  ^  h,  which  is  consistent  with  Equation  (14).  Further¬ 


more,  the  number  of  triangles  in  To  that  have  an  edge  on  any  one  side  of  [0,  If  is  given 
by 


It  follows  that  there  are  ( Ne  +  l)2  nodes  and  2N^  triangles  in  Tq. 


4.1  Construction 

The  Cartesian  coordinates  of  the  nodes  that  constitute  the  vertices  of  the  triangles  in  To 
are  generated  by  CoarseNodes  (Algorithm  4)  and  stored  in  AC  The  order  of  the  nodes  in 
J\f  is  shown  in  Figure  2. 
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Figure  2:  A  coarse  triangulation  showing  the  assignment  of  indices  to  nodes  and  triangles, 
and  the  grid  spacing  5. 


Algorithm  4:  CoarseNodes 

Input:  Ne  £  N 
Output:  A f 

1  5  -  N-1 

2  N «—  (vM[l},...,VM[{Ne  +  l)2]) 

3  k  <—  1 

4  foreach  i  €  (0, 1, ,  Are)  do 

5  foreach  j  £  (0, 1, ... ,  Ne)  do 

6  N[k]  <-  6  (j,  i) 

7  k  <—  k  +  1 

8  return  J\f 


The  vertices  of  the  triangles  in  To  are  stored  in  7^v  as  references  to  their  unique  entries 
in  J\f .  References  to  the  vertices  of  the  zth  triangle  are  stored  in  7^v  [i] :  the  order  of  the 
Tq  [i]  agrees  with  the  assignment  of  indices  to  triangles  shown  in  Figure  2.  Since  L{t)  =  0 
for  all  t  £  To,  the  order  of  the  three  references  in  each  Tq  [i]  corresponds  to  an  anticlockwise 
arrangement  of  the  vertices  such  that  the  first  and  last  vertices  form  the  bisection  edge; 
see  Algorithm  1.  The  data  structure  Tq  is  constructed  by  CoarseTriangleVertices  as 
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per  Algorithm  5.  Note  that  only  contains  references  to  vertices  and  therefore  Tq  is 
initialised  using  the  null  data  structure  given  by  Vj-  =  ((0,  0,  0)  , . . . 


Algorithm  5:  CoarseTriangleVertices 


Input:  Ne  £  N 
Output:  Tq 

i  n  <— 

Ag+1 

2  Tq  * 

-(T°[l],... 

,T°[2IVe2]) 

3  k 

<- 

1 

4  foreach  i  S  (0,1, 

. . .,  Ne  -  1)  do 

5 

foreach  j  G  (0, 1, . . . ,  Ne  —  1)  do 

6 

T0V[M]- 

j  +  1  +  n(i  +  1) 

7 

VM- 

j  +  1  +  ni 

8 

VM]- 

j  +  2  +  ni 

9 

k  <—  k  +  1 

10 

T0V[M]- 

j  +  2  +  ni 

11 

Vik,  2]  - 

j  +  2  +  n  (i  +  1) 

12 

VM]  - 

■  j  +  1  +  n(i  +  1) 

13 

k  <—  k  +  1 

14  return  Tq 

The  neighbours  of  the  ith  triangle  in  To  are  stored  in  '7qN  [?i]  as  references  to  their 
unique  triangle  indices;  the  assignment  of  indices  to  triangles  is  shown  in  Figure  2. 
The  order  of  the  three  references  in  each  '7qN  [i]  coincides  with  the  attributes  of  T  for 
storing  triangle  neighbours;  see  Section  2.1.  The  data  structure  7^N  is  constructed  by 
CoarseTriangleNeighbours  in  accordance  with  Algorithm  6,  where  the  index  (3  <  0  de¬ 
notes  a  reference  to  a  null  boundary  triangle.  Observe  from  Figure  2  that  triangles  with 
an  odd  index  may  have  a  boundary  edge  on  the  left  and  bottom  segments  of  the  boundary 
of  To,  whereas  triangles  with  an  even  index  may  have  a  boundary  edge  on  the  right  and 
top  segments  of  the  boundary  of  To;  only  the  first  and  last  triangles  have  two  boundary 
edges. 
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Algorithm  6:  CoarseTriangleNeighbours 

Input:  Ne  G  N 
Output:  7^n 

1  n  <—  2 IVg  //  number  of  triangles  in  To 

2  s  <—  2Ae  //  number  of  triangles  in  a  strip 

3  Z  < —  1  //  left  boundary  indicator 

4  r  <—  1  //  right  boundary  indicator 

5  A;  «—  n  —  s  +  1 

6  T0N^  (U0[l],...,U0[n]) 

7  T0N[1]  <-  (2,P,P) 

8  T0N[n]  <-  (n  -  1 ,/?,/?) 

9  foreach  *  G  (2, 3, . . . ,  n  —  1)  do 

10  if  i  is  odd  then 

n  if  *  =  Is  +  1  then 

12  l  < —  l  1 

13  7^N[z]  ^  (i  +  l,/3,  z  -  s  +  1) 

14  else  if  i  >  s  then 

15  I  2oN [*]<—(*+ M  —  M  —  s  +  1) 

16  else 

17  L  roNN  <-  (*+  M  -  hP) 

18  else 

19  if  i  =  rs  then 

20  r  <—  r  +  1 

21  7^N[i]  <—  (i  -  l,/3,  i  +  s  -  1) 

22  else  if  i  <  k  then 

23  |  7^N[z]  <— (i  —  l,z  +  l,z  +  s  —  1) 

24  else 

25  |_  r0N[z]  <-  (i  -  l,i  +  l,P) 

26  return  7^,N 
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4.2  Mesh  Initialisation 

The  results  of  the  previous  sections  can  now  be  employed  to  initialise  T.  First  To  is 
constructed  in  accordance  with  Algorithms  4  to  6  and  stored  in  the  data  structures  J\f  and 
T ;  the  variables  IVX  and  Nt  are  also  stored.  Equations  (2)  and  (4)  are  then  used  to  estimate 
the  number  of  triangles  in  the  final  refined  T,  and  the  amount  of  null  data  required  to  store 
T  is  appended  to  M  and  T .  This  procedure  is  executed  by  InitialiseMesh  (Algorithm  7), 
which  has  O(Nt)  running  time. 


Algorithm  7:  InitialiseMesh 
Input:  S  >  0  and  h(  ■ )  as  per  Equation  (16) 

Result :  T  is  initialised  and  stored  in  J\[  and  T ;  Ax  and  Nt  are  also  stored 

1  Ne  <  [1/5] 

2  IVX  <-  (JVe  +  l)2 

3  Nt  «-  2IVg 

4  N  <—  CoarseNodes(IVe) 

5  Tq  <—  CoarseTriangleVertices(Ae) 

6  7()N  <—  CoarseTriangleNeighbours(IVe) 

7  T^(Vr[l],...,Vr[Nt\) 

8  foreach  i  €  (1,2,...,  Nt)  do 

9  |_  T[i]  e-  (T0V [i,  1] ,  T0V [i,  2] ,  Tq  [z,  3] ,  T0N [z,  1] ,  T0N [z,  2] ,  T0N [z,  3] ,  1 , 0) 

10  6  <-  1/Ne 

U  12eb1(S’Jl(i)) 

12  Ni  e-  [Nl/2\ 

13  N  <-  (AT[1] ,  N[2] , . . . ,  AT[IVX] ,  PM[1] , . . . ,  TV [iv£  -  Nx]) 

14  T^(T[l],T[2},...,T[Nt\,VT[l},...,Vr[N^-Nt}) 


Algorithms  5,  6  and  the  initialisation  of  T  on  line  9  of  Algorithm  7  establish  the 
attributes  of  T  discussed  in  Section  2.1.  Also  note  that  the  target  local  triangulation 
diameter  function,  given  by  Equation  (12),  is  modified  to  take  a  triangle  index  as  an 
argument  for  use  on  line  11  of  Algorithm  7: 

h{i)=  min  h{N[j]),  (16) 

jeT[i,  1:3] 


for  i  €  {1,2, . .  .,Nt}. 
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5  Mesh  Refinement 

This  section  concludes  the  implementation  of  Maubach’s  method.  Algorithms  for  simulta¬ 
neously  bisecting  compatibly  divisible  triangles  and  updating  their  neighbours  are  exam¬ 
ined;  these  are  the  most  important  algorithms  of  the  implementation  as  their  behaviour 
is  critically  linked  to  Maubach’s  method  and  the  data  structures  M  and  T .  The  top- 
level  mesh  refinement  algorithm  is  also  presented,  together  with  a  demonstration  of  its 
computational  performance. 


5.1  Bisecting  Triangles 

BisectTriangle  bisects  a  given  triangle  ti  €  T  and  creates  the  two  descendants.  A  high- 
level  description  of  BisectTriangle  is  given  in  Algorithm  1.  The  objective  of  this  section 
is  to  provide  a  detailed  implementation  of  an  algorithm  that  simultaneously  bisects  two 
compatibly  divisible  triangles  in  accordance  with  Maubach’s  method. 

It  will  be  shown  in  Section  7  that  point  location  can  be  achieved  using  a  constant 
running-time  algorithm  without  storing  descendants.  Therefore,  information  regarding 
tj.  is  not  retained  after  its  bisection  and  the  index  i  is  reused  to  reference  one  of  the 
descendants  of  U.  The  other  descendant  of  U  is  given  the  index  i  =  Nf  +  1,  where  Nt 
represents  the  state  of  Nt  immediately  prior  to  bisecting  ti. 

Maubach’s  method  stipulates  that  the  vertices  of  descendants  be  ordered  to  enable  the 
descendants’  bisection  edges  to  be  subsequently  identified  without  any  computations.  If 
U  is  the  triangle  to  be  bisected,  this  order  depends  on  the  value  of  k(tt).  where 

k{ti)  =  2  —  ( L(ti )  mod  2)  €  {1,2}. 

Let  V(ti)  =  (a,  b ,  c).  If  k{ti)  =  1,  then  a  and  b  reference  the  nodes  that  form  the  bisection 
edge,  otherwise  a  and  c  reference  the  nodes  that  form  the  bisection  edge. 

Consider  bisecting  ti  for  the  case  k{ti)  =  1  and  O(U)  =  0,  which  corresponds  to  a 
clockwise  arrangement  of  the  vertices.  Then  by  Maubach’s  method,  the  vertices  of  the 
descendants  L  and  tj  will  satisfy 


V(U)  =  (a,d,c),  (17) 

V(ti)  =  (b,d,c),  (18) 

where  the  new  node  with  index  d  is  stored  in  A f[Nx  +  1],  that  is,  d  =  1VX  +  1.  The 
configuration  of  the  descendants  U  and  tj  is  shown  in  Figure  3  for  the  case  being  considered. 
Observe  from  Figure  3  and  Equations  (17)  and  (18)  that  the  descendants  satisfy  0(ti)  =  0 
and  0{tj)  =  1. 

This  example  reveals  the  dependence  of  the  configuration  of  the  descendants  of  ti  on 
k(ti)  and  O(L).  The  configuration  of  descendants  must  be  known  to  update  triangle 
neighbours  without  any  computations,  which  is  discussed  in  Section  5.2.  Since  k(U)  and 
0(ti)  can  each  take  two  values,  there  are  four  configurations  of  the  descendants.  The 
remaining  three  configurations  can  be  extrapolated  from  Figure  3  by  ordering  the  vertices 
according  to  Maubach’s  method. 
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a  a 


Figure  3:  The  configuration  of  the  descendants  ofti  for  the  case  k(ti )  =  1  and  Oftf)  =  0. 
The  references  to  the  descendants  are  i  and  i,  and  a,  b,  c  and  d  denote  references  to 
vertices.  Note  that  Nx  and  Nt  represent  their  states  immediately  prior  to  bisecting  ti. 


Let  i  be  the  index  of  a  triangle  to  be  bisected  and  %  the  index  of  its  bisection  neigh¬ 
bour.  If  tib  is  not  a  null  boundary  triangle  and  ti  and  tib  are  compatibly  divisible,  then 
BisectTriangles(i)  (Algorithm  8)  simultaneously  bisects  ti  and  t,;fj ,  creating  their  de¬ 
scendants  in  accordance  with  Maubach’s  method,  and  updates  their  neighbours. 

The  bisection  of  a  triangle  with  its  bisection  edge  on  the  boundary  of  D  is  undertaken  by 
BisectBoundaryTriangle  (Algorithm  19  in  Appendix  B).  A  compatibly  divisible  triangle 
with  a  non-bisection  edge  on  the  boundary  of  D  is  bisected  using  BisectTriangles. 

Note  that  for  the  algorithms  in  Section  5  and  Appendix  B,  the  symbols  Nx  and  Nt 
appearing  in  the  “Input”  and  “Result”  fields  represent  the  states  of  Nx  and  Nt  immediately 
prior  to  calling  the  routines. 
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Algorithm  8:  BisectTriangles 

Input:  i  €  {1,2, . . . ,  Nt}  such  that  T[i,  4]  >  0  and  T[i,  8]  =  T[ib,  8] 

Result:  T[i\  and  T[ib\  are  updated;  A /”[AX  +  1],  T[Nt  +  1]  and  T[Nt  +  2]  are 
created 


1  CheckMemory( ) 

2  %<-  T[i,  4] 

3  V  <-  T[i,  1  :  3] 

4  Vb  *—  T[ib,  1  :  3] 

5  A  x  <—  Ax  +  1 

6  k(i )  <—  2  —  (T[i,  8]  mod  2) 

7  Af[Ax]  <-  0.5  (A^[U[1]]  +  M[V[k  +  1]]) 


8  Nt  < —  Nf  +  2 

9  O  <—  (T[z,  7]  +  1)  mod  2 

10  Ob  <—  (' T[ib,  7]  +  1)  mod  2 

11  T[i,  8]  «-  T[i,  8]  +  1 

12  T[i6, 8]  <—  T[ib,  8]  +  1 


13  if  k(i )  =  1  then 


14 

T[i,  1], 

-V[l] 

15 

T[i,2]  < 

-Ax 

16 

T[i,  3]  < 

-U[3] 

17 

T[%,  1] 

-  Vb[l] 

18 

T\ib,  2] 

Ax 

19 

T[ib,  3] 

<-  H[3] 

20 

T[Nt- 

l]^(U[2],Ax,U[3],0,0,0,0,r[i,8]) 

21 

T[Nt\  <- 

-  (Vb{2\ ,  Ax,  Vb[3] ,  0, 0, 0,  Ob,  T\ib,  8]) 

22  else 


23 

24 

25 

26 

27 

28 

29 

30 


T[i,l]^V[l] 
m  2]+-V[2] 

T\i ,  3]  <-  IVx 

r[i&,  i]  <-  vb[i\ 

T[ib,2]  <-  Vb[2] 

T[ib,  3]  <—  Ax 

T[Nt  -  1]  <-  (U[2],  U[3],  Ax,  0, 0, 0,  T[i,  7],  T[i,  8]) 
T[Nt]  <-  (I4[2],  Vb[3],  Ax,  0, 0, 0,  T[ib,  7],T[ib,  8]) 


31  switch  (k(i),  T[i,  7],  T[ib,  7])  do 

32  case  (1,  0,  0)  or  (2, 1, 1) 

33  |  UpdateNeighboursl(i) 

34  case  (1, 1,  0)  or  (2,  0, 1) 

35  |  UpdateNeighbours2(i) 

36  case  (1,  0, 1)  or  (2, 1,  0) 

37  |  UpdateNeighbours3(i) 

38  otherwise 

39  UpdateNeighbours4(i) 
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5.2  Updating  Triangle  Neighbours 

It  was  shown  in  Section  5.1  that  the  configuration  of  the  descendants  of  ti  can  be  de¬ 
termined  from  the  values  of  kitf)  and  O(fj)  prior  to  bisection.  This  will  now  be  used  in 
conjunction  with  the  attributes  of  T  to  update  the  references  to  neighbours  of  compatibly 
divisible  triangles  post  bisection,  without  any  computations. 

Consider  bisecting  compatibly  divisible  triangles  ti  and  tib  for  the  case  k(U)  =  1, 
0(ti)  =  0  and  0(Ub )  =  0.  The  attributes  of  T  stipulate  that  ib  and  i  be  the  first  references 
to  the  neighbours  of  L  and  tib,  respectively.  Let  the  second  and  third  references  to  the 
neighbours  of  tt  and  tib  satisfy 


N(U)  =  ( ib,a,b ), 


(19) 


N(tib)  =  ( i,c,d ). 


(20) 


The  configuration  of  ti  and  Ub  and  their  neighbours  prior  to  bisection  is  shown  on  the 
left  of  Figure  4;  this  follows  from  Equations  (19)  and  (20)  and  the  attributes  of  T.  The 
configuration  of  the  descendants  of  ti  and  tib  and  their  neighbours  is  shown  on  the  right  of 
Figure  4;  this  follows  from  Maubach’s  method  as  per  Figure  3.  Observe  that  ti  and  tb  are  no 
longer  neighbours  after  bisection;  likewise,  tib  and  are  no  longer  neighbours.  Therefore 
N(tb)  and  N(td)  must  also  be  updated  to  reflect  these  changes;  this  is  performed  by  call¬ 
ing  ReplaceTriangleNeighbour  (Algorithm  20  in  Appendix  B).  UpdateNeighboursl(i), 
which  is  called  by  BisectTriangles  (Algorithm  8),  updates  the  neighbours  of  the  descen¬ 
dants  of  ti  and  tib  for  the  case  being  considered.  An  implementation  of  UpdateNeighboursl 
is  presented  in  Algorithm  9. 


d 


> 


d 


Figure  f:  The  configuration  of  triangle  neighbours  before  and  after  the  bisection  of  com¬ 
patibly  divisible  triangles  ti  and  tib  for  the  case  k{tf)  =  1,  Oftf)  =  0  and  0(tib )  =  0.  Here 
a,  b,  c  and  d  represent  triangle  neighbours. 


This  example  reveals  the  dependence  of  the  configuration  of  neighbours  on  k(ti),  Ofti ) 
and  0(tib).  Since  each  of  these  can  take  two  values,  there  are  potentially  eight  config¬ 
urations  of  neighbours.  However  it  can  be  shown  that  only  four  of  these  configurations 
are  independent.  The  remaining  three  independent  configurations  can  be  extrapolated 
from  Figures  3  and  4  by  employing  Maubach’s  method  and  the  attributes  of  T.  Imple¬ 
mentations  of  UpdateNeighbours  for  these  cases  can  be  found  in  Algorithms  21  to  23  in 
Appendix  B. 
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Algorithm  9:  UpdateNeighboursl 

Input:  i  €  {1,2, . . .  ,Nt  —  2}  such  that  (fc(i),  T[i,  7],  T[4,  7])  =  (1, 0,  0)  or  (2, 1, 1) 
Result:  T[i,  4  :  6]  and  T[p,,  4  :  6]  are  updated;  T[Nt  —  1,4  :  6]  and  T[IVt,4  :  6]  are 
created 

1  i+-  Nt-  1 

2  ib<~  Nt 

3  ib<-  T[i,  4] 

4  A  <—  T[z,  5  :  6] 

5  Nb<-  T[ib ,  5  :  6] 

6  T[i,4]  4-  IV[1] 

7  T[i,  5]  <—  i 

8  T[i,  6]  4-  ib 

9  T[i6,4]4-JVb[l] 

10  T[ib,  5]  <—  ib 

11  T[ib,  6]  <—  i 

12  T[i,  4]  e-  IV[2] 

13  T[i,  5]  4—  ib 

14  T[i,  6]  4—  i 

15  T[ib,  4]  4-  Nb[ 2] 

16  T[i;,,  5]  4—  i 

17  T[i6, 6]  4-  ib 

18  ReplaceTriangleNeighbour  ( IV [2] ,  i,  i ) 

19  ReplaceTriangleNeighbour (iVj,[2] ,  ib,  ib) 


Consider  bisecting  a  triangle  U  with  its  bisection  edge  on  the  boundary  of  fl.  The 
configuration  of  the  descendants  of  L  and  their  neighbours  only  depends  on  k(ti)  and 
0(ti),  as  ti  does  not  have  a  neighbour  on  its  bisection  edge  prior  to  bisection.  The  neigh¬ 
bours  of  the  descendants  of  L  are  updated  by  calling  UpdateBoundaryNeighboursl  or 
UpdateBoundaryNeighbours2,  depending  on  the  values  of  k(tj)  and  O(L);  implementa¬ 
tions  of  these  routines  can  be  found  in  Appendix  B. 
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5.3  Refinement  Algorithm 

An  updated  version  of  Ref  ineTriangle  that  incorporates  the  algorithms  of  the  previous 
sections  is  presented  in  Algorithm  10.  Most  of  the  challenges  of  implementing  Maubach’s 
method  originate  in  the  data  structures  for  storing  T  and  the  routines  for  updating  triangle 
neighbours.  This  is  reflected  by  the  similarity  of  Algorithm  10  and  the  more  abstract 
version  of  Ref  ineTriangle  appearing  in  Algorithm  2. 

Ref  ineTriangle(i)  calls  itself  repeatedly,  starting  with  ti,  until  a  compatibly  divisible 
triangle  tc  is  found;  see  Algorithm  10.  This  generates  a  sequence  of  triangles.  A  stack 
(last-in-first-out  queue)  Sb  enables  the  triangles  in  this  sequence  to  be  bisected  in  reverse 
order  from  tc  to  ti.  Therefore  compatibility  is  preserved  because  only  compatibly  divisible 
triangles  are  bisected;  this  is  performed  by  Algorithm  11.  RefineMesh  (Algorithm  12) 
initialises  Sb  then  calls  Ref  ineTriangle  repeatedly  until  the  local  refinement  condition 
(on  line  5  of  Algorithm  12)  returns  False  for  all  t  £  T,  and  then  removes  the  excess  null 
data  from  J\f  and  T . 

Finally,  a  refined  mesh  on  [0,  l]2  is  generated  using  Maubach’s  method  by  specifying 
the  maximum  grid  spacing  5  and  the  target  local  triangulation  diameter  h(  ■ ),  and  calling 
InitialiseMesh(<5,  h(  ■ ))  followed  by  Ref  ineMesh(c),  h(  ■ )).  Extending  this  implementa¬ 
tion  to  allow  for  non-square  domains  and  domains  with  holes  is  considered  in  Appendix  A. 


Algorithm  10:  Ref  ineTriangle 


Input:  i  <E  { 1,2, . . . ,  Nt} 
Result:  ti  is  compatibly  refined 


l  ib  <-  T[i,  4] 


2 

3 

4 

5 

6 

7 

8 
9 

10 


if  ib  <  0  then 

BisectBoundaryTriangle(i) 
BisectTrianglesInStack( ) 
else  if  T[i,  8]  =  T[ib,8\  then 
BisectTriangles(i) 
BisectTrianglesInStack( ) 
else 

push  i  onto  Sb 
Ref ineTriangle(ife) 


Algorithm  11:  BisectTrianglesInStack 
Result:  triangles  in  stack  Sb  are  compatibly  bisected 

1  while  Sb  is  not  empty  do 

2  j  <—  pop  the  top  element  from  Sb 

3  BisectTriangles(j) 
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Algorithm  12:  RefineMesh 

Input:  5  >  0  and  h(  ■ )  as  per  Equation  (16) 
Result:  the  final  refined  T  is  stored  in  J\f  and  T 

1  initialise  S 5 

2  <5^1/  \l/I\  //  see  Equation  (15) 

3  i  <—  1 

4  while  i  ^  Nt  do 

5  while  £b(8,T[i,  8])  >  h(i)  do 

6  |_  Ref  ineTriangle(i) 

7  i  <—  i  +  1 

8  Af  J\f[l  :  Nx] 

9  T  <-  T[1  :  Nt] 


5.4  Computational  Performance 

It  was  noted  in  Section  1.1  that  Maubach’s  method  has  O(Ni)  running  time.  It  will  now  be 
demonstrated  that  0(Nt )  running  time  can  be  attained  in  practice  when  InitialiseMesh 
and  RefineMesh  are  used  for  mesh  generation. 

The  mesh  shown  in  Figure  1  was  used  for  testing  the  computational  performance  of 
InitialiseMesh  and  RefineMesh.  The  maximum  grid  spacing  was  fixed  at  8  =  0.05  and 
the  target  local  triangulation  diameter  was  chosen  to  be 

_  r(i,  <0.4 

I  1,  otherwise, 

where  0  <  d  <  \/2  5.  Values  of  d  were  chosen  that  resulted  in  the  final  refined  trian¬ 
gulations  containing  between  9,560  and  2,113,944  triangles,  with  a  maximum  13  levels  of 
refinement.  Average  CPU  times  and  a  least-squares  line  of  best  fit  were  calculated  after 
running  InitialiseMesh  and  RefineMesh  10  times  for  each  value  of  d7 

Figure  5  shows  the  scaled  average  CPU  times  f  =  t/tq^q  versus  Nt  =  Nt/ 9560,  where 
T9560  is  the  average  CPU  time  to  generate  the  refined  triangulation  with  9,560  triangles. 
The  O(Nt)  running  time  of  InitialiseMesh  and  RefineMesh  for  the  triangulations  being 
considered  is  evident  in  Figure  5,  where  the  line  of  best  fit  has  a  gradient  of  approximately 
1.01. 


'95%  confidence  intervals  were  also  calculated,  but  they  were  too  small  to  be  visible  on  the  plot. 
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Figure  5:  The  scaled  average  running  time  of  the  implementation  of  Maubach ’s  method 
documented  in  this  report,  versus  the  scaled  number  of  triangles;  a  line  of  best  fit  is  shown 
with  the  data  points. 


6  Adjacency  List 


Finding  the  neighbours  of  a  given  node  is  a  common  operation  on  a  mesh.  The  neighbours 
of  a  node  x  are  the  nodes  that  are  adjacent  to  x;  two  nodes  are  adjacent  if  they  form  an 
edge  of  a  triangle.8  An  adjacency  list  A  is  an  ordered  collection  of  unordered  lists,  where 
A[i\  is  the  set  of  all  node  indices  that  reference  the  neighbours  of  the  node  x,-. 


The  adjacency  list  can  be  constructed  during  mesh  generation  in  an  analogous  manner 
to  the  initialisation  and  maintenance  of  the  mesh  data  for  triangle  neighbours;  see  Sec¬ 
tions  4.1  and  5.2.  Alternatively,  the  adjacency  list  can  be  constructed  from  an  existing 
mesh;  this  approach  is  considered  here. 


AdjacencyList  (Algorithm  13)  constructs  A  by  finding  all  the  edges  in  T.  Edges 
are  constructed  from  triangle  vertex  references  by  visiting  each  t  £  T  once;  this  procedure 
duplicates  references  to  node  neighbours  as  each  edge  is  shared  by  two  triangles.9  Therefore 
each  list  in  A  must  have  capacity  for  16  indices,  since  a  node  can  have  at  most  eight 
neighbours  in  a  mesh  generated  using  Maubach’s  method.  After  finding  all  the  edges  in 
T,  the  duplicated  indices  and  any  excess  padding  are  deleted  from  each  list  in  A. 


sThe  neighbours  of  x  can  also  be  given  by  the  set  yK(x)  \  x,  where  JV (x)  is  given  by  Equation  (9). 

9  Boundary  edges  are  an  exception. 
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Algorithm  13:  AdjacencyList 


Input:  Nx  and  Tv 
Output:  A 


1  Pv  * 

-((2, 3), (1,3), (1,2)) 

2  pA  * 

—  (!)••■)  //  \PA\=NX 

3  A  <- 

-  ((0, 0, 0,  0, 0,  0,  0,  0,  0,  0,  0,  0,  0, 0,0,0),.. 

. )  //  \A\  =  Nx 

4 

5 

6 

7 

8 
9 

10 


foreach  i  G  (l,  2, . . . ,  |TV|)  do 
foreach  j  G  (1,2,3)  do 

<-  p.aM 

A[n,np]  <-  Tv[i,pv[j,  1]] 
A[n,  np  +  1]  <-  Tv [i,Pv[j,  2]] 
P.aM  np  +  2 


n  foreach  jg  (1,2,...,  Nx)  do 

12  A[i]  <—  delete  duplicates  from  A[i] 

13  A[i]  <—  delete  zero  from  A[i] 

14  return  A 


AdjacencyList  takes  Nx  and  Tv  as  inputs,  where  Tx  is  a  data  structure  contain¬ 
ing  references  to  the  triangle  vertices  of  T;  Tv[i]  =  T[i,  1  :  3]  for  i  G  {1,2,...,  -ZV*}. 
References  to  the  vertices  that  form  unordered  edges  of  all  t  G  T  are  obtained  from 
Pv  =  ((2,  3),  (1,3),  (1,2))  in  Algorithm  13,  that  is,  (Tv[i,  2],  Tv[i,  3]),  (Tv  [i,  1],  Tv[i,  3]) 
and  (Tv[i,  1],  Tv  [i,  2])  are  references  to  the  unordered  edges  of  The  list  in  Algo¬ 
rithm  13  enables  node  neighbours  to  be  added  to  the  lists  in  A  via  an  in-place  change:  the 
next  neighbour  of  node  i  to  be  found  is  placed  in  position  pj\\i]  of  A[i\.  Most  programming 
languages  provide  functions  for  efficiently  deleting  duplicates  from  a  list,  and  hence  the 
details  of  line  12  of  Algorithm  13  are  omitted. 

The  implementation  of  AdjacencyList  given  by  Algorithm  13  has  the  following  prop¬ 
erties: 

•  A  is  constructed  in  0(NX)  running  time,  as  the  two  loop  bodies  each  have  0(1) 
running  time 

•  0(NX )  memory  is  required  to  store  A 

•  references  to  node  neighbours  are  obtained  from  A  in  fast  0(1)  running  time. 


7  Point  Location 

Given  a  query  point  y  G  f2,  find  at  G  T  that  geometrically  contains  y.  This  is  known  as  the 
point  location  problem.  Point  location  is  a  fundamental  operation  on  a  mesh  and  therefore 
has  many  applications  [Miicke,  Saias  &  Zhu  1999,  Soukal,  Maikova  &  Kolingerova  2012]. 
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One  prominent  application  of  point  location  is  interpolation.  Let  g(x)  be  a  quantity 
that  is  known  at  each  x  6  'U(T),  and  suppose  there  is  a  requirement  to  calculate  q(y)  for 
some  y  6  Si.  After  using  a  point  location  method  to  find  a  i  6  T  that  contains  y,  ^(y ) 
can  be  interpolated  from  the  vertices  Xj  6  t  using,  for  example,  barycentric  co-ordinates: 

Ci  +  C2  +  C3  =  i, 

y  =  Cixi  +  C2X2  +  C3X3, 

q(  y)  =  Ci(?(xi)  +  C2<?(x2)  +  C3  ?(x3). 

Point  location  is  accomplished  in  two  steps: 

1.  select  a  t  6  T  to  begin  the  search 

2.  traverse  T  until  a  t  £  T  containing  the  query  point  y  is  found. 

A  judicious  selection  of  the  initial  triangle  is  critical  for  an  efficient  implementation  of  point 
location.  A  dedicated  data  structure  can  be  used  to  perform  the  second  step;  a  binary 
tree  is  an  obvious  choice  for  a  mesh  generated  using  Maubach’s  method.  Alternatively, 
the  existing  triangle  neighbour  data  can  be  employed  to  “walk”  from  triangle  to  triangle 
until  a  triangle  containing  y  is  found.  The  two  steps  of  point  location  are  examined  in  the 
following  two  sections. 


7.1  Initial  Triangle  Selection 

A  technique  that  subdivides  D  into  “buckets”  using  a  uniform  grid  can  be  used  to  select 
the  initial  triangle  by  identifying  a  bucket  containing  y  [Asano  et  al.  1985].  The  simplicity 
of  this  subdivision  enables  a  bucket  containing  y  to  be  found  in  constant  running  time, 
however  additional  memory  is  required  to  store  the  buckets. 

An  alternative  to  bucketing  that  does  not  require  additional  memory  begins  by  ran¬ 
domly  choosing  a  subset  of  points  from  Y{T).  A  point  from  this  subset  that  minimises 
the  distance  to  y  is  then  selected  to  be  the  starting  point  of  a  straight-line  walk  to  y.  The 
resulting  point  location  algorithm  has  an  0(N x'  )  expected  running  time  [Miicke,  Saias 
&  Zhu  1999]. 

For  a  mesh  generated  using  Maubach’s  method,  the  initial  triangle  can  be  selected 
in  constant  running  time  without  requiring  additional  memory.  This  follows  from  three 
factors: 

•  Maubach’s  method  uses  bisection  refinement  and  therefore  all  descendants  of  a  coarse 
triangle  are  geometrically  contained  within  that  triangle 

•  data  regarding  the  history  of  descendants  is  not  stored  and  the  index  of  a  parent 
triangle  is  reused  to  reference  one  of  its  descendants 

•  Algorithm  7  generates  a  uniform  triangulated  grid  Tq  to  initialise  Maubach’s  method. 
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Hence  if  i  is  the  index  of  a  coarse  triangle  in  To  that  contains  y,  then  the  vertices  of  T  G  T 
are  at  most  a  distance  of  h  from  y. 

Fast  initial  triangle  selection  for  point  location  is  the  key  reason  Algorithm  7  was  chosen 
to  initialise  Maubach’s  method  in  preference  to  simply  using  two  compatible  triangles. 
Algorithm  7  enables  fast  initial  triangle  selection  by  creating  virtual  buckets.  Therefore 
initial  triangle  selection  can  be  achieved  by  finding  a  coarse  triangle  in  To  containing  y; 
this  is  performed  by  InitialTriangle  in  Algorithm  14. 


Algorithm  14:  InitialTriangle 

Input:  Ne  G  N  and  y  €  [0,  l]2 

Output:  i  G  {1, 

1  5  <-  1/Ne 

2  (x,y)  <-  y 

3  x'  <—  max  {1,  [ Ne  x\ } 

4  y'  i  max { 1 ,  \Ney]} 

5  *  «-  2x'  -  1  +  2Ne(y'  -  1) 

6  if  y  >  5(x '  +  y'  —  1)  —  x  then 

7  |  return  i  +  1 

8  else 

9  return  i 


InitialTriangle  finds  a  coarse  triangle  in  To  that  contains  y  in  constant  running 
time  by  exploiting  the  simple  geometry  of  uniform  grids  and  the  numbering  of  triangles 
shown  in  Figure  2.  First  a  square  containing  y  is  ascertained  by  counting  the  triangles 
that  would  be  traversed  if  the  search  were  to  begin  in  U  G  To,  followed  by  U  G  Tj,  and 
so  on:  on  line  5  of  Algorithm  14,  2 Ne(y'  —  1)  is  the  number  of  coarse  triangles  in  each 
traversed  row,  and  2x'  —  1  is  the  minimum  number  of  coarse  triangles  traversed  from  the 
left  of  the  grid  to  a  square  containing  y.  Then  a  coarse  triangle  containing  y  will  be  above 
or  below  the  line  connecting  the  upper  left  and  lower  right  corners  of  a  square  containing 
y;  this  test  is  performed  on  line  6  of  Algorithm  14. 


7.2  Triangle  Containing  the  Query  Point 

After  employing  InitialTriangle  to  initialise  point  location,  T  is  traversed  until  a  t  G  T 
is  found  that  contains  the  query  point  y.  This  section  is  devoted  to  the  examination  of 
a  walking  method  for  traversing  T  that  uses  the  orientation  test  shown  in  Figure  6  to 
determine  the  next  step  in  the  walk.  The  PointLocation  routine  (Algorithm  15)  utilizes 
this  walking  method  to  find  a  t  containing  y. 

The  objective  of  a  walking  method  is  to  move  closer  to  y  with  each  step  in  the  walk. 
Let  t  be  the  current  triangle  being  tested  for  containing  y.  PointLocation  effectively 
bisects  D  along  each  edge  of  t  and  then  steps  into  the  first  neighbour  of  t  that  belongs  to 
the  subdomain  containing  y.  The  orientation  test  determines  the  subdomain  containing 
y:  if  t  does  not  contain  y,  then  y  will  be  to  the  right  of  one  of  the  edges  of  i;  see  Figure  6. 
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a 


J 


Figure  6:  The  orientation  test:  does  d  lie  on,  to  the  left  of,  or  to  the  right  of  the  line 
defined  by  the  triangle  edge  (a,  b )  ?  The  arrows  indicate  the  directions  of  the  vectors  used 
to  perform  the  test.  The  test  is  applied  to  the  other  edges  in  the  same  manner. 


If  y  is  not  to  the  right  of  any  edge  of  t,  then  t  must  contain  y  and  PointLocation  returns 
the  index  of  t  and  terminates. 

The  orientation  test  is  performed  by  establishing  the  sign  of  the  determinant  of  a 
matrix  with  entries  constructed  from  the  edge  and  point  being  tested  [Soukal,  Maikova  & 
Kolingerova  2012].  For  the  example  shown  in  Figure  6,  this  is  equivalent  to  determining 
the  sign  of  the  third  component  of  the  cross  product  (d  —  b)  x  (a  —  b). 10  The  orientation 
test  is  performed  on  line  14  of  Algorithm  15. 

The  orientation  test  is  known  to  produce  incorrect  results  when  y  and  the  edge  being 
tested  are  approximately  colinear  [Shewchuk  1997a].  A  detailed  examination  of  this  issue 
is  beyond  the  scope  of  this  report;  however  the  impact  of  an  erroneous  orientation  test 
is  examined  here.  Let  the  edge  e  of  t  be  approximately  colinear  with  y,  and  assume  the 
orientation  test  has  just  produced  an  incorrect  result.  Consider  these  three  possibilities: 

1.  a  wrong  neighbour  of  t  is  stepped  into,  or 

2.  subsequent  steps  cycle  between  neighbouring  triangles  that  do  not  contain  y,  or 

3.  subsequent  steps  oscillate  between  a  triangle  and  its  neighbour  that  does  contain  y. 

The  first  problem  is  inconsequential  as  the  outcome  will  simply  be  a  somewhat  longer  walk. 
The  second  problem  is  solved  by  testing  the  edges  of  t  in  a  random  order  for  each  step 
of  the  walk  [Devillers,  Pion  &  Teillaud  2002],  This  is  successful  because  cycles  can  only 
occur  if  PointLocation  continues  to  erroneously  step  into  the  neighbour  of  t  that  shares 
e:  since  y  cannot  be  colinear  with  two  edges  of  t  simultaneously,  y  must  be  to  the  right 
of  one  edge  that  is  not  e,  and  therefore  testing  this  edge  before  e  will  halt  the  cycle.  The 
third  problem  is  solved  by  preventing  PointLocation  from  returning  to  the  triangle  that 
was  previously  tested,  which  will  also  improve  the  efficiency  of  PointLocation  [Devillers, 
Pion  &  Teillaud  2002].  However,  this  may  result  in  PointLocation  returning  a  triangle  t! 
that  does  not  contain  y.  Even  so,  in  this  case  y  must  be  approximately  on  the  edge  of  t! 

10The  points  a ,  b  and  d  are  considered  here  to  be  3-dimensional  vectors  in  the  xy  plane. 
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shared  with  its  neighbour  that  does  contain  y,  and  the  resulting  numerical  error  may  be 
tolerable  for  the  application  employing  point  location. 

PointLocation  requires  the  vectors  used  to  perform  the  orientation  test  to  have  an 
anticlockwise  orientation  about  the  triangle  being  tested,  as  per  Figure  6.  The  data 
structure  £[i]  is  used  to  ensure  the  edges  of  t%  form  vectors  with  the  desired  orientation, 
in  an  efficient  manner.11  The  edges  of  all  L  £  T  are  represented  implicitly  by  £[i],  where 


'((1,2),  (3,1),  (2,  3)) 
((2,1), (3, 2), (1,3)) 
((3,1), (2, 3), (1,2)) 
k  ((1,3),  (2,1),  (3,  2)) 


k{i)  =  1  and  0(i)  =  0 
k(i)  =  1  and  0(i)  =  1 
k(i)  =  2  and  0(i)  =  0 
k(i)  =  2  and  0(i)  =  1, 


for  *£{1,2,...,  Ntj.  £[i]  has  the  following  attributes: 


•  £[*,  j]  represents  the  jth  edge  of  U 

•  £[i,  1]  represents  the  bisection  edge  of  L  with  the  remaining  edges  arranged  anti¬ 
clockwise  about  ti 

•  £[i,j]  for  j  =  1,2,3  (in  order)  are  the  positions  of  T\i,  1  :  3]  containing  references  to 
nodes  that  enable  vectors  to  be  formed  that  rotate  anticlockwise  about  tt. 


For  example,  consider  the  case  k(i)  =  1  and  O(i)  =  0.  Then  Af[T[i,  1]]  —  Af[T[i,2]], 
A/"[T[*,3]]  —  J\f[T[i,  1]]  and  A/"[T[*,2]]  —  jV[T[i,  3]]  (in  order)  are  vectors  that  rotate  anti¬ 
clockwise  about  ti . 

An  implementation  of  PointLocation  is  presented  in  Algorithm  15.  PointLocation 
calls  InitialTriangle  and  returns  the  index  of  a  triangle  containing  the  query  point  y. 
Algorithm  15  has  constant  running  time  and  requires  insignificant  additional  memory  to 
store  £.  Recall  that  if  i  =  InitialTriangle(lVe, y),  then  the  vertices  of  L  G  T  are  at 
most  a  distance  of  h  from  y.  Equivalently,  a  triangle  containing  y  will  be  at  most  0{ 2Ml) 
steps  away  from  L ,  where  Ml  =  maxJ-er12r..jivt}  L(tj).  Hence  the  number  of  steps  taken 
by  PointLocation  is  limited  to  2Ml .  The  orientation  test  is  performed  on  line  14  of 
Algorithm  15. 

Note  that  Algorithm  15  is  only  valid  for  H  =  [0,  l]2.  Extending  Algorithm  15  to  allow 
for  non-square  domains  and  domains  with  holes  is  considered  in  Appendix  A. 


1£  can  also  be  used  with  the  boundary  marker  (3  to  efficiently  determine  the  references  to  nodes  on  the 
boundary  of  Q.  This  is  particularly  useful  for  a  general  Q,  which  is  considered  in  Appendix  A. 
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Algorithm  15:  PointLocation 


Input:  JVe£N  and  y  G  D 

Output:  i  €  {1, 


prev 


o 


o 


"-max 

B  «-  True 

i  <—  InitialTriangle(lVe,  y) 

6  t  <-  (Af[T[i,l}},  N[T[i,2]],  JV[T[i,3]]) 

7  while  B  and  n  ^  nmax  do 
B  <—  False 

9  E  <-  £[i] 

10  P  <—  random  permutation  of  {1,2,3} 

n  for  each  j  e  P  do 

12  (n,r2)  e- y  -  t[E[j,2]] 

13  (ei,e2)  <-  t[E[j,  1]]  -  t[E[j,2]\ 

//  do  the  orientation  test,  but  don’t  go  back  or  step  out  of  D 
if  T[i,j  +  3]  /  *prev  and  T[i,  j  +  3]  >  0  and  sign(uie2  —  r2ei)  >  0  then 

B  <-  True 


14 

15 

16 

17 

18 

19 

20 


6 prev  T  1 

i  J  +  3] 

t  (Af[T[i,l}},  Af[T[i,2}},  Ar[T[i,3}}) 

break 


n  <—  n  +  1 


21  return  i 


8  Conclusion 


This  report  provides  a  comprehensive  implementation  of  the  unstructured  mesh  genera¬ 
tion  method  of  Maubach  [1995],  focussing  on  the  case  of  a  two-dimensional  mesh  on  a 
square  domain.  An  extension  to  this  implementation  that  enables  mesh  generation  on 
two-dimensional  non-square  domains  and  domains  with  holes  is  presented  in  Appendix  A. 

The  implementation  has  the  following  desirable  features.  Mesh  data  structures  were 
chosen  to  enable  local  bisection  refinement  to  occur  in  constant  running  time  and  with 
minimal  computations.  A  local  refinement  condition  was  derived  that  guarantees  the  local 
triangulation  diameter  of  the  refined  mesh  will  obey  a  specified  bound;  hence  the  numer¬ 
ical  error  of  computations  on  the  mesh  can  be  controlled  while  restricting  the  number  of 
triangles.  The  mesh  refinement  algorithm  was  tested  and  shown  to  achieve  the  anticipated 
linear  running  time  with  respect  to  the  number  of  triangles  in  the  refined  mesh.  Employ¬ 
ing  a  uniform  triangulated  grid  to  initialise  Maubach’s  method  creates  virtual  buckets. 
It  follows  that,  for  a  mesh  generated  using  Maubach’s  method,  point  location  can  be 
accomplished  in  constant  running  time  without  requiring  additional  memory. 
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Appendix  A  Domains  other  than  [0, 1] 

The  implementation  of  Maubach’s  method  documented  in  the  body  of  this  report  is  only 
valid  for  D  =  [0, 1] 2 .  This  implementation  is  now  extended  to  allow  for  non-square  do¬ 
mains  and  domains  with  holes,  by  covering  the  general  domain  D  with  a  square  domain, 
initialising  Maubach’s  method  using  this  square  domain,  and  then  removing  the  triangles 
not  in  D.  Maubach’s  method  can  then  be  used  to  refine  the  triangles  in  Q. 

First  consider  fl  =  [a,b]2.  Let  Da 0  =  b  —  a  and  yc  =  (2/1 , 2/2)  such  that 

^  =  [yi  -  Doo/2,  yi  +  -Doo/2]  x  [y2  -  D0 c/2,  y2  +  D^/2]  .  (Al) 

Let  u:  [0,  l]2  — >  [a,  b]2  where 

u(y)  =  Doo  (y  -  (0.5,  0.5))  +  yc.  (A2) 

The  inverse  of  u  is  u~1 :  [a,  b]2  — >  [0,  l]2  where 

^_1(y)  =  vv—  (y  -  yc)  +  (0.5, 0.5).  (A3) 

In  this  case  Ne  and  5  are  given  by: 

Ne  =  \°f],  (A4) 

s=^-  <A5> 

The  following  changes  must  be  made  to  accommodate  the  case  D  =  [a,  b]2: 

1.  use  the  definition  of  Ne  given  by  Equation  (A4)  in  all  algorithms 

2.  use  the  definition  of  5  given  by  Equation  (A5)  in  all  algorithms 

3.  immediately  after  line  4  of  Algorithm  7,  replace  J\f  with  the  result  of  applying  the 
function  u  to  each  element  of  J\f 

4.  replace  y  with  u_1(y)  on  line  5  of  Algorithm  15. 

Now  consider  the  general  case  where  D  is  non-square,  possibly  with  holes.  Assume 
there  exists  a  routine  DomainQ(y)  that  returns  True  if  y  G  12,  or  False  otherwise,  and 
determine  a  square  domain  as  per  Equations  (Al)  to  (A5)  that  covers  D,  that  is,  if  y  €  D 
then  y  £  [a,6]2.  To  accommodate  the  general  case,  make  the  above  changes  for  a  square 
domain  and  modify  InitialiseMesh  (Algorithm  7)  as  follows: 

5.  immediately  after  line  6  of  Algorithm  7,  call  FilterNodes  (Algorithm  16) 

6.  then  call  FilterVertices  (Algorithm  17) 

7.  then  call  FilterNeighbours  (Algorithm  18). 

Finally, 

8.  replace  line  5  of  Algorithm  15  with  i  <—  Mt  [lnitialTriangle(fV4,  u_1(y))] , 
where  Mt  is  constructed  by  Algorithm  17. 
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Algorithm  16:  FilterNodes 
Input:  AT  generated  by  Algorithm  4 

Result:  nodes  not  in  D  are  removed  from  AT,  Afx  is  constructed  and  IVX  is  updated 

1  3  <-  1 

2  A4X  <—(1,2,...,  IVX) 

3  AT  <—  (7VI1],  ■  •  •  i  'Pm[Nx\) 

4  n  <—  Nx 

5  IVX  < —  0 

6  foreach  i  e  (1,  2, . . . ,  n)  do 

7  if  DomainQ(AA[?'])  then 

8  Nx  <—  Nx  +  1 

9  A/’fA'x]  <-  Af[i] 

10  A4x[i]=j 

11  j  <-  j  +  1 

12  else 

13  |_  Afx[i]  =  0 

14  Af  <—  A/"[l  :  Ax] 


Algorithm  17:  FilterVertices 
Input:  7^v  generated  by  Algorithm  5 

Result:  references  to  nodes  not  in  D  are  removed  from  Tq  ,  A4 1  is  constructed  and 
Nt  is  updated 

1  1 

3  Mt  <-  (1,2,...,  Nf) 

4  n  <-  A{ 

5  lVt<-0 

//  triangles  that  do  not  have  all  three  vertices  in  fl  are  removed 

6  foreach  i  e  (1,  2, . . . ,  n)  do 

7  V^(MX  [T0V [i,  1]] ,  Mx  [T0V [i,  2]] ,  Mx  [T0V [i,  3]] ) 

8  if  U[l]  >  0  and  U[2]  >  0  and  U[3]  >  0  then 

9  N^+-Nt  +  1 

10  VlNt}^-  V 

11  A4t[i]=j 

12  j  «-  j  +  1 

13  else 

14  |_  AU[«]  =  /3 

15  T0V  e-  tJ[1  :  At] 
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Algorithm  18:  FilterNeighbours 

Input:  7q^  generated  by  Algorithm  6 

Result:  references  to  triangles  not  in  Q  are  removed  from  7^,N 

1  k  <—  0 

2  n  <-  |T0N| 

3  N  <-  (0,0,0) 

4  ],..., U°[n]) 

//  n  ^  Nt  when  FilterNeighbours  is  called 

5  foreach  i  E  (1,  2, . . . ,  n)  do 

6  if  AU[i]  >  0  then 

7  foreach  j  E  (1,2,3)  do 

8  if  7)(n  [i,  j]  >  0  then 

9  |  N[j}^Mt[T"[i,j}\ 

10  else 

n  [  L  "W-r0N[i,j] 

12  k  <—  k  +  1 

13  |_  T^[k]  =N 

//  now  k  =  Nt 

14  T0n  e-  7*[1  :  L] 
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Appendix  B  Other  Cases  for  Bisecting  Triangles 
and  Updating  their  Neighbours 


Algorithm  19:  BisectBoundaryTriangle 

Input:  i  G  {1,2, . . . ,  Nt}  such  that  T[i,  4]  <  0 

Result:  T[i]  is  updated;  A 7[AX  +  1]  and  T[Nt  +  1]  are  created 

1  CheckMemory( ) 

2  V  <-  T[i,  1  :  3] 

3  A'x  Ax  +  1 

4  k(i )  2  —  (T[z,  8]  mod  2) 

5  A r[Nx\  <-  0.5  (Af[U[l]]  +  A f[V[k  +  1]]) 

6  Nt  <-  Nt  +  1 

7  O  < —  (T[i,  7]  +  1)  mod  2 
s  T[i,  8]  <-  T[i,  8]  +  1 


9  if  k  =  1  then 


10 

T[i,  1]  <■ 

-U[l] 

11 

T[i,2]e 

-Ax 

12 

T[i,  3]  e 

-U[3] 

13 

T[Nt]  e 

-(U[2],ATX,U[3],0,0,0,0,T[A8]) 

14  else 

15 

T[i,  1]  <■ 

-U[l] 

16 

T[i,  2]  e 

~V[2] 

17 

T[i,  3]  <■ 

-  Ax 

18 

_  T[Nt]  e 

-  (^[2]  T[3] ,  Arx,  0, 0, 0,  T[i,  7] ,  T[i,  8]) 

19  switch  (k{i),T[i,  7])  do 


20 

21 

22 

23 


case  (1, 0)  or  (2, 1) 

j  UpdateBoundaryNeighboursl(i) 
otherwise 

I-  UpdateBoundaryNeighbours2(i) 


Algorithm  20:  ReplaceTriangleNeighbour 


Input:  i,  j,k  <E  {1,2, ...  ,Nt} 

Result:  ti,  which  was  a  neighbour  of  tj,  is  updated  to  be  a  neighbour  of  tk 


l  if  i  >  0  then 


if  T[i,  4]  =  j  then 
|  T[i,  4]  <-  k 
else  if  T\i,  51  =  j  then 
I  T[i,  5]  «-  k 
else 

L  6]  <—  k 
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Algorithm  21:  UpdateNeighbours2 

Input:  i  €  {1,2, . . .  ,Nt  —  2}  such  that  (fc(i),  T[i,  7],  T[4,  7])  =  (1, 1,  0)  or  (2,  0, 1) 
Result:  T[i,  4  :  6]  and  T[ib,4  :  6]  are  updated;  T[Nt  —  1,4  :  6]  and  T[Ntl  4  :  6]  are 
created 

1  i  <-  Nt  -  1 

2  ib<~  Nt 

3  ib  <-  T[z,4] 

4  IV  <—  T[i,  5  :  6] 

5  Ab  <-  T[ib, 5  :  6] 

6  T[i,4]  «-  A[2] 

7  T[i,  5]  <-  ib 

8  T[i,  6]  <—  i 

9  T[i6,4]  Ab[l] 

10  T[ift,  5]  «-  ib 

11  T[ib,  6]  «-  i 

12  T[i,  4]  «-  A[l] 

13  T[i,  5]  <—  i 

14  T[i,  6]  <-  ib 

is  T|z6,4]<-  JVfe[2] 

16  T[4, 5]  <—  i 

17  T[ift,  6]  «-  ib 

18  ReplaceTriangleNeighbour  ( A[l] ,  i,  i ) 

19  ReplaceTriangleNeighbour ( JVj,[2] ,  ib ,  ib) 
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Algorithm  22:  UpdateNeighbours3 

Input:  i  €  {1,2, . . . ,  Nt  —  2}  such  that  (k(i),T[i,  7],  T[ib,  7])  =  (1,  0, 1)  or  (2, 1,  0) 
Result:  T[i,  4  :  6]  and  T[ib,  4  :  6]  are  updated;  T[Nt  —  1,4  :  6]  and  T[A^,4  :  6]  are 
created 

1  i<-  Nt-1 

2  ib<~  Nt 

3  ib  <-  T[z,4] 

4  A  <—  T[i,  5  :  6] 

5  Ab  <-  T[zb,  5  :  6] 

6  T[i,  4]  <-  A[l] 

7  T[i,  5]  <—  i 

8  T[i,  6]  <-  ib 

9  T[i6,4]  <-  Ab[2] 

10  5]  <—  i 

11  T[ib,  6]  «-  ib 

12  T[i,  4]  «-  A[2] 

13  T[i,  5]  <—  ib 

14  T[i,  6]  <—  i 

is  T\ibA\  <-  Nb[l] 

16  T[ib,  5]  <-  ip 

17  T[ib,  6]  «-  i 

18  ReplaceTriangleNeighbour  ( A [2] ,  i,  i ) 

19  ReplaceTriangleNeighbour (A^[l] ,  ib,  ib) 
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Algorithm  23:  UpdateNeighbours4 

Input:  i  €  {1,2, . . .  ,Nt  —  2}  such  that  (fe(i),  T[i,  7],  T[4,  7])  =  (1, 1, 1)  or  (2,  0,  0) 
Result:  T[i,  4  :  6]  and  T[ib,  4  :  6]  are  updated;  T[A7  —  1,4  :  6]  and  T[Nt,  4  :  6]  are 
created 

1  i  <-  AT*  -  1 

2  ib<~  Nt 

3  ib  <-  T[z,4] 

4  IV  <—  T[i,  5  :  6] 

5  Ab  <-  T[i6, 5  :  6] 

6  T[i,4]  <-  A[2] 

7  T[i,5]  <—  ib 

8  T[i,  6]  <—  i 

9  T[ib,4]  Ab[2] 

10  5]  <—  i 

11  T[ib,  6]  «-  ih 

12  T[i,  4]  «-  A[l] 

13  T[i,  5]  <—  i 

14  T[i,6]  <-  ib 

is  T[i6,4]  <-  Nb[l] 

16  T[ib,  5]  <-  ib 

17  6]  «-  i 

18  ReplaceTriangleNeighbour  ( A[l] ,  i,  i) 

19  ReplaceTriangleNeighbour ( Nb[l] ,  ib,  ib ) 
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Algorithm  24:  UpdateBoundaryNeighboursl 

Input:  i  G  {1,2, . . .  ,Nt  —  1}  such  that  (k(i),T[i,  7])  =  (1,0)  or  (2, 1) 
Result:  T[i,  4  :  6]  is  updated;  T[Nt,4  :  6]  is  created 

i  i^Nt 
2/?<-T[i,4] 

3  N  <—  T[i,  5  :  6] 

4  T[i,  4]  <-  IV[1] 

5  T[i,  5]  <—  i 

6  T[i,  6]  <—  f3 

7  T\l ,  4]  «-  IV[2] 

8  T[i,  5]  <—  /3 

9  T[i,  6]  <—  i 

10  ReplaceTriangleNeighbour  ( iV[2] ,  i,  i) 


Algorithm  25:  UpdateBoundaryNeighbours2 

Input:  i  G  {1,2, . . .  ,Nt  —  1}  such  that  (fe(i),  T[i,  7])  =  (1, 1)  or  (2,  0) 
Result:  T[i,  4  :  6]  is  updated;  T[Nt,4  :  6]  is  created 

i  i*-Nt 
2f3^T[i,4] 

3  N  <—  T[i,  5  :  6] 

4  T[i,  4]  <-  A[2] 

5  T[i,  5]  <—  f3 

6  T[i,  6]  <—  i 

7  T[i,  4]  <-  A[l] 

8  T[i,  5]  <—  i 

9  T[i,  6 ]  *—  (3 

10  ReplaceTriangleNeighbour (A[l] ,  i,  i) 
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